# -*- coding: utf-8 -*-
"""
Created on Tue May  7 13:25:21 2019

@author: Solydo
"""
import numpy as np

def softmax(X,Y,lamuda=0.01,alpha=0.1,e=1E-5):
    # softmax 回归，梯度下降法
    # 超参数： lamuda-规则化系数，alpha-学习率，e-收敛条件
    # 返回参数theta，设有N个类别，格式为[t1,t2,...tN],ti为n+1维列向量，最后一个数值为截距b
    labels=list(set(Y))  #所有可能的类标记
    N=len(labels)        #类别数目
    m,n=np.shape(X)       #m为样本数，n为特征数
    X1=np.c_[X,np.ones([m,1])]  #X1在X基础上增加一个全1列,维度变为 m×(n+1)
    Y1=np.zeros([m,N])          #Y1为扩展后的Y，维度m×C
    for c in range(N):
        Y1[Y==labels[c],c]=1 
    theta=np.zeros([n+1,N])    #初始化theta参数，维度(n+1)×C
    dis=1  #梯度大小
    loop=1
    while dis>e and loop<5E3:
        p=np.exp(np.dot(X1,theta))
        p=p/sum(p.T).reshape(-1,1)  #维度m×C,计算每个样本属于每个类别的概率
        d=np.dot(X1.T,(Y1-p))/(-m)+lamuda*theta   #梯度
        dis=np.sqrt(sum(d.reshape(-1)**2))      #梯度的模长
        theta=theta-alpha*d
        loop+=1
    return theta

def AnalysisData(X,DataType=None):
    # 分析数据
    # DataType 返回每个属性的数据类型，若为字符串类型，说明是离散型，并统计各个可能取值
    # X1 返回数字化的数据，将离散型数值转变为相应取值的编号
    m=len(X)     #样本数
    n=len(X[0])  #属性数
    if DataType==None: #如果未从外部输入DataType,则需自动创建
        DataType=[]
        for p in range(n):
            DataType.append([type(X[0][p])])
            if DataType[p][0]==str:
                for q in range(m):
                    if X[q][p] not in DataType[p]:
                        DataType[p].append(X[q][p])
    X1=[]   #数字化的数据
    for p in range(m):
        temp=[]
        for q in range(n):
            if DataType[q][0]==str:
                temp.append(DataType[q].index(X[p][q]))
            else:
                temp.append(X[p][q])
        X1.append(temp)
    return DataType,X1

def most(Y):
    # 统计样本集中的多数类
    if len(Y)==0:
        return None
    Yvalue=[]  #Y取值
    Ynum=[]    #取值数目
    for y in Y:
        if y not in Yvalue:
            Yvalue.append(y)
            Ynum.append(1)
        else:
            Ynum[Yvalue.index(y)]+=1
    return Yvalue[Ynum.index(max(Ynum))]

def Index(Y,rule='InfoGain'):
    # 计算信息熵或者基尼系数
    if len(Y)==0:
        if rule=='InfoGain':
            return 0
        elif rule=='Gini':
            return 1
    Yvalue=[]  #Y取值
    Ynum=[]    #取值数目
    for y in Y:
        if y not in Yvalue:
            Yvalue.append(y)
            Ynum.append(1)
        else:
            Ynum[Yvalue.index(y)]+=1
    pk=np.array(Ynum)/sum(Ynum)
    if rule=='InfoGain':
        return sum(-pk*np.log2(pk))
    elif rule=='Gini':
        return 1-sum(pk**2)  

def score(X,Y,considerfeature,DataType,rule='InfoGain'):
    # 这里considerfeature为单个考察属性的索引号
    # 仅适用于信息熵和基尼系数两种选择规则，哪种规则由参数rule确定
    # 由于计算信息熵和基尼值的过程多有相似，因此统一编程
    # 返回信息熵/基尼值
    # 若为连续属性，还要返回划分点
    f=considerfeature
    if DataType[f][0]==str:  #若该特征取值为离散型
        Gain=Index(Y) #信息增益初始值
        Gini=0
        for v in range(1,len(DataType[f])):
            Yv=Y[X[:,f]==v]
            Gain-=Index(Yv)*len(Yv)/len(Y)
            Gini+=Index(Yv,'Gini')*len(Yv)/len(Y)
        Gain=[Gain]
        Gini=[Gini]
    else:   #若该特征取值为连续型
        Gain=[]
        Gini=[]
        values=np.array(list(set(X[:,f])))  #利用集合去除重复取值
        if len(values)==1:  #若当前子集X中该特征取值相同
            Gain=[0,values]
            Gini=[Index(Y,'Gini'),values]
        else:
            t=(np.sort(values)[:-1]+np.sort(values)[1:])/2  #候选划分点
            for tt in t:
                Y1=Y[X[:,f]<=tt]; Y2=Y[X[:,f]>tt]
                Gain.append(Index(Y)-Index(Y1)*len(Y1)/len(Y)-Index(Y2)*len(Y2)/len(Y))
                Gini.append(Index(Y1,'Gini')*len(Y1)/len(Y)+Index(Y2,'Gini')*len(Y2)/len(Y))
            Gain=[max(Gain),t[Gain.index(max(Gain))]]
            Gini=[min(Gini),t[Gini.index(min(Gini))]]
    if rule=='InfoGain':
        return Gain
    if rule=='Gini':
        return Gini

def DividePoint(X,Y,lamuda):
    # 输入X仅含单个属性
    # 针对对率回归作为选择规则时，对于连续属性，确定其划分点
    # 若仅含两个类别，直接将对率回归得到的分界面作为划分点
    # 若多于两类，比如N类，将得到N组theta值，两两一对可以确定出一个分界面
    # 多个分界面，考察其训练误差，选择训练误差最小的那个

    theta=softmax(X,Y,lamuda=0.01)
    N=theta.shape[-1]       #类别数
    if N==2:                #若仅含两个类别
        dt=theta[:,0]-theta[:,1]
        return -dt[1]/dt[0]
    else:
        t=[]       #划分点
        right=[]   #正确分类数
        for i in range(N):
            for j in range(i+1,N):
                dt=theta[:,i]-theta[:,j]
                t.append(-dt[1]/dt[0])       #划分点
                #print(X,Y,t)
                Y1=Y[X[:,0]<=t[-1]]
                Y2=Y[X[:,0]>t[-1]]
                right.append(sum(Y1==most(Y1))+sum(Y2==most(Y2)))  #训练正确分类数
        return t[right.index(max(right))]
            
def OneHot(X,DataType,order):
    # one-hot独热编码，将原数据X变为独热编码的数据X1
    if order==None:
        order=[True]*len(DataType)
    num=[len(i)-1 for i in DataType]  #每个属性的取值数
    expand=[num[i]*(not order[i])+order[i] for i in range(len(num))] #每个属性需要扩展的数目
    X1=[]
    for i in range(len(num)):
        if expand[i]==1:
            X1=np.r_[X1,X[:,i]]
        else:
            for n in range(expand[i]):
                X1=np.r_[X1,(X[:,i]==(n+1))]
    return X1.reshape(sum(expand),-1).T

def TransConsider(considerfeature,DataType,order):
    #======== 输入 ========
    # 当前考察属性(属性索引号) considerfeature
    # 各个属性的数据信息 DataType
    # 各个属性的有序性 order
    #======== 输出 ========
    # 对无序属性的数据进行one-hot编码之后
    # 所考察属性新的索引号 newcf
    # 新的索引号所对应的原索引号 IndexBack
    num=[len(i)-1 for i in DataType]  #每个属性的取值数
    expand=[num[i]*(not order[i])+(order[i]) for i in range(len(num))] #每个属性需要扩展的数目
    newcf=[]
    IndexBack=[]
    for cf in considerfeature:
        IndexBack+=[cf]*expand[cf]
        newcf+=list(np.arange(expand[cf])+sum(expand[:cf]))
    return newcf,IndexBack

def choose(X,Y,FeatureName,considerfeature,DataType,rule='InfoGain',lamuda=0,equalchoose=False,order=None):
    # 选择最佳划分属性，基于三种选择规则：信息增益、基尼系数、对率回归
    scores=[] #下面还有values用于存储各个属性的
              #信息增益、基尼系数或者对率回归的权重平均序值
              #scores除了存储这些指标，还有划分点信息(对于连续取值属性)
    if rule=='LogisticReg':
        global scale,bias
        X1=(OneHot(X,DataType,order)-np.array(bias))/np.array(scale)  #归一化
        newcf,IndexBack=TransConsider(considerfeature,DataType,order) #所考察属性的索引号需要转换
        theta=softmax(X1[:,newcf],Y,lamuda=0.01)    #softmax回归拟合到的theta参数
        values=np.mean(abs(theta[:-1,:]),axis=1)    #各个属性的平均权重值
        values=list(values)
        bestvalue=max(values)                       #最大的权重值
        index=values.index(bestvalue)               #最大权重值对应的属性索引号
        f=IndexBack[index]                          #上面索引号是对于one-hot编码后的，需要转换回去
        if DataType[f][0]==str:                     #若为离散型属性，只需返回属性索引号
            return [f]
        else:
            t=DividePoint(X[:,[f]],Y,lamuda)        #若为连续型属性，还要返回划分点 
            return [f,t]
    if (rule=='InfoGain')|(rule=='Gini'):
        for f in considerfeature:
            scores.append(score(X,Y,f,DataType,rule))
        values=[v[0] for v in scores]
    if rule=='InfoGain':
        bestvalue=max(values)
    if rule=='Gini':
        bestvalue=min(values)

    #当多个属性值对应的信息增益/基尼系数相等时，提示选择其中一个 
    if values.count(bestvalue)!=1 and equalchoose:
        index=np.arange(len(values))[values==bestvalue]
        print('现有%d个属性所得信息增益/基尼系数/对率回归权重序值相等,请从中选择一个：'%len(index))
        k=1
        for i in index:
            print('%d:%s'%(k,FeatureName[considerfeature[i]]),end=', ')
            k+=1
        index=index[eval(input())-1]
    else:
        index=values.index(bestvalue)
    
    if len(scores[index])>1:  #返回最佳属性编号，若为连续取值，还要返回划分点
        return [considerfeature[index],scores[index][1]]
    else:
        return [considerfeature[index]]

def precut(Xt,Yt,Xv,Yv,BestFeature,DataType,equalcut=True):
    # 预剪枝子函数
    right0=sum(Yv==most(Yt)) #划分前分类正确数,(多数类分类正确)
    right1=0                 #划分后分类正确数
    bf=BestFeature[0]
    if DataType[bf][0]==str:  #若划分属性为离散属性
        for k in range(1,len(DataType[bf])):
            Ytk=Yt[Xt[:,bf]==k]
            Yvk=Yv[Xv[:,bf]==k]
            if len(Yvk)!=0:
                if len(Ytk)!=0:
                    right1+=sum(Yvk==most(Ytk))
                else:
                    right1+=sum(Yvk==most(Yt))
    else:   #若划分属性为连续属性
        t=BestFeature[1]
        Yt1=Yt[Xt[:,bf]<=t]
        Yv1=Yv[Xv[:,bf]<=t]
        if len(Yv1)!=0:
            if len(Yt1)!=0:
                right1+=sum(Yv1==most(Yt1))
            else:
                right1+=sum(Yv1==most(Yt))
        Yt2=Yt[Xt[:,bf]>t]
        Yv2=Yv[Xv[:,bf]>t]
        if len(Yv2)!=0:
            if len(Yt2)!=0:
                right1+=sum(Yv2==most(Yt2))
            else:
                right1+=sum(Yv2==most(Yt))
    cut=(right0>right1)|((right0==right1) & equalcut)
    #若划分后正确数低于划分前，则不进行划分，亦即预剪枝
    #若划分前后正确数相等，要是同时参数equalcut=True,也要剪枝
    #=====================================
    #以下语句调试时用，用于显示具体计算到的划分前后的正确数
    '''
    print('当前划分特征的序号:',BestFeature,',划分前:',right0,',划分后:',right1)
    if cut:
        print('预剪枝, 不划分')
    else:
        print('不预剪枝, 要划分')
    '''
    #=====================================
    return cut

def aftercut(tree,Xt,Yt,Xv,Yv,FeatureName,LabelName,DataType,equalcut=False):
    # 后剪枝子函数,递归方法实现
    #--------------计算剪枝后正确分类数-------------
    right1=sum(Yv==most(Yt))  #剪枝后的正确分类数(将树变为叶结点，类别为训练集多数类)
    #--------------计算剪枝前正确分类数-------------
    right0=0
    if type(tree)!=dict:  #若为单个叶结点
        return tree,right1
    NodeText=list(tree.keys())[0]              #NodeText为结点文本
    f=FeatureName.index(NodeText.strip('=?'))  #f为所考察属性的索引号
    discrete=(DataType[f][0]==str)             #当前考察属性是否为离散值
    subs=tree[NodeText]                        #subs为分支内容，字典类型
    for key in subs.keys():  #具体每个分支的分支信息
        if discrete:
            k=DataType[f].index(key)   #将分支取值'key'转换为序号第k个
            Xtk=Xt[Xt[:,f]==k,:]
            Ytk=Yt[Xt[:,f]==k]
            Xvk=Xv[Xv[:,f]==k,:]
            Yvk=Yv[Xv[:,f]==k]
        else:
            Xtk=eval('Xt[Xt[:,f]'+key+',:]')
            Ytk=eval('Yt[Xt[:,f]'+key+']')
            Xvk=eval('Xv[Xv[:,f]'+key+',:]')
            Yvk=eval('Yv[Xv[:,f]'+key+']')            
        if type(subs[key])==dict: #如果该分支含有子树,则对子树进行后剪枝，并返回剪枝结果和正确分类数
            CuttedTree,right=aftercut(subs[key],Xtk,Ytk,Xvk,Yvk,FeatureName,LabelName,DataType,equalcut)
            tree[NodeText][key]=CuttedTree #修改该分支内容
            right0+=right
        else:#如果该分支为叶结点
            # ↓该叶结点所属类号
            y=list(LabelName.keys())[list(LabelName.values()).index(subs[key])]
            right0+=sum(Yvk==y)
    #剪后正确数增加了,或者不变且要剪
    cut=(right1>right0)|((right1==right0) and equalcut)
    #=====================================
    #下面语句用于调试时显示详细的剪枝详细
    '''
    print('考察属性:',NodeText,',剪枝前:',right0,'剪枝后:',right1)
    if cut:   
        print('后剪枝:剪了')
    else:
        print('后剪枝:不剪')
    '''
    #=====================================
    if cut:
        return LabelName[most(Yt)],right1
    else:
        return tree,right0

def predict(tree,X,FeatureName,LabelName):
    # 决策树对样本集X进行预测，采用递归方法实现
    # 输入X为列表类型，数据可以为数值或者字符串
    # 输出Y为列表类型，代表预测结果
    m=len(X)
    Y=[0]*m
    if type(tree)!=dict: #当决策树为单结点时
        return [list(LabelName.keys())[list(LabelName.values()).index(tree)]]*m
    NodeText=list(tree.keys())[0]              #NodeText为结点文本
    f=FeatureName.index(NodeText.strip('=?'))  #f为所考察属性的索引号
    discrete=(type(X[0][f])==str)              #当前考察属性是否为离散值
    subs=tree[NodeText]                        #subs为分支内容，字典类型
    for key in subs.keys():  #具体每个分支的分支信息
        index=[]             #属于该分支的样本索引号
        if discrete:
            for i in range(m):
                if X[i][f]==key:
                    index.append(i)
        else:
            for i in range(m):
                if eval('X[i][f]'+key):
                    index.append(i)
        if len(index)==0:
            continue
        if type(subs[key])==dict: #如果该分支含有子树
            Ysub=predict(subs[key],[X[i] for i in index],FeatureName,LabelName)
            for i in range(len(index)):
                Y[index[i]]=Ysub[i]
        else: #如果该分支为叶结点
            # ↓该叶结点所属类号
            y=list(LabelName.keys())[list(LabelName.values()).index(subs[key])]
            for i in range(len(index)):
                Y[index[i]]=y
    return Y

scale=[]  #***  全局变量，用于对数据进行归一化操作
bias=[]   #***  X'=(X-bias)/scale

def CreatTree(Xt,Yt,Xv,Yv,FeatureName,LabelName,considerfeature='all',DataType=None,rule='InfoGain',cut=None,equalcut=True,deep=0,normal=None,lamuda=0.01,equalchoose=False,order=None):
    #================================================================
    #                      递归方法训练决策树
    # *****输入变量*****
    # Xt,Xv---------训练集和剪枝集的X---二层列表数据类型，列表内容可以是字符串或者数值
    # Yt,Yv---------训练集和剪枝集的X---单层列表数据类型，列表内容为数值
    # FeatureName---特征(属性)名称。----单层列表数据类型，列表内容为字符串                
    # LabelName-----类标记名称。--------单层字典类型，key为数值型，对应于Yt/Yv中所有可能数值，value为字符串
    # consi***ture--当前考察特征--------单层列表类型，内容为所考察特征所对应的索引号，默认为‘all’，意为所有特征
    # DataType------数据类型------------二层列表类型，格式为比如：[[type('str'),'a','b'],[type(0.0)]]
    #                                      意为有两个特征，第一个特征数据为字符串类型，可能取值为'a'和'b',
    #                                      第二个特征为数值型。
    # rule----------决策树划分规则------字符串类型，可能三种取值:'InfoGain'-信息熵,'Gini'-基尼,'LogisticReg'-对率回归
    # cut-----------剪枝类型------------字符串类型，可能三种取值：'pre'-预剪枝,'after'-后剪枝，None-不剪枝
    # equalcut------精度相等时是否剪枝---逻辑类型:True/False
    # deep----------当前决策树深度-------整数类型，初始值为0，(实际上相当递归深度)
    # normal--------归一化方式-----------字符串类型，三种取值：'min-max','z-score',None
    # lamuda--------规则化参数-----------数值类型，仅用于对率回归
    # equalchoose---是否手动选择属性-----逻辑类型:True/False，用于多个属性的信息熵/基尼值/权重值相等时
    # order---------各个属性的有序性-----单层列表类型，内容和True/False,表明各个属性的有序性
    # *****输出变量*****
    # tree----------决策树--------------多层嵌套字典类型
    #================================================================  
    #print(deep,considerfeature)
    #-------------初始化--------------
    if considerfeature=='all':   #初始时的'all'表示考虑所有特征
        considerfeature=list(range(len(FeatureName)))
    if deep==0:                  #初始时需要分析数据，将X数字化，并得到每个X的数据类型
        if order==None:
            order=[True]*len(FeatureName)
        mt=len(Yt)
        mv=len(Yv)
        DataType,X=AnalysisData(Xt+Xv,DataType)
        X=np.array(X)
        Xt=X[:mt,:]
        Xv=X[mt:,:]
        Yt=np.array(Yt)
        Yv=np.array(Yv)
        
        global scale,bias     #确定对率回归时的归一化操作参数 (X-bias)/scale ->X'
        X1=OneHot(X,DataType,order) #独热编码
        if normal==None:      #不做归一化
           bias=[0]*X1.shape[1]
           scale=[1]*X1.shape[1]
        if normal=='min-max': #min-max归一化
           bias=list(np.min(X1,axis=0))
           scale=list(np.max(X1,axis=0)-np.min(X1,axis=0))
        if normal=='z-score': #z-score归一化
           bias=list(np.mean(X1,axis=0))
           scale=list(np.std(X1,axis=0))
    #------ --考察几种特殊情况---------
    if (Yt==Yt[0]).all():      #若所有样本属于同一类
        return LabelName[Yt[0]]
    if len(considerfeature)==0:  #若待考察特征为空
        return LabelName[most(Yt)]
    if (Xt[:,considerfeature]==Xt[0,considerfeature]).all():  #若所有样本取值相同
        return LabelName[most(Yt)]
    #--------最优划分属性和预剪枝-------
    BestFeature=choose(Xt,Yt,FeatureName,considerfeature,DataType,rule,lamuda,equalchoose,order)  #选择最佳划分属性
    bf=BestFeature[0]
    #print(deep,BestFeature)
    if cut=='pre' and precut(Xt,Yt,Xv,Yv,BestFeature,DataType,equalcut):
        return LabelName[most(Yt)]
    #------------分支-----------------
    if len(BestFeature)==1:  #若为离散型属性
        nt=FeatureName[bf]+'=?'  #nt意为nodetext
        tree={nt:{}}
        considerfeature.remove(bf)
        for k in range(1,len(DataType[bf])):
            Xtk=Xt[Xt[:,bf]==k,:]
            Ytk=Yt[Xt[:,bf]==k]
            Xvk=Xv[Xv[:,bf]==k,:]
            Yvk=Yv[Xv[:,bf]==k]
            if len(Ytk)==0:  #若该分支的训练集为空
                tree[nt][DataType[bf][k]]=LabelName[most(Yt)]
            else:
                tree[nt][DataType[bf][k]]=CreatTree(Xtk,Ytk,Xvk,Yvk,FeatureName,LabelName,considerfeature.copy(),DataType,rule,cut,equalcut,deep+1,normal,lamuda,equalchoose,order)
    else: #若为连续型属性
        nt=FeatureName[bf]  #nt意为nodetext
        tree={nt:{}}
        t=BestFeature[1]    #划分点
        t=float('%.4f'%t)   # 这一句是考虑到这样一种情况：
                            # 比方t=0.999935,划分子集时，x=1将被划分到'x>t'
                            # 然而在决策树中存储的信息为t=1.0,见下面语句tree[nt]['<=%.3f'%t]=...
                            # 在后剪枝或者预测新样本时，x=1将被划分到'x<=t'分支
                            # 将会出现前后不一致的情况，甚至于出错
        #----------'<=t'分支-----------
        Xt1=Xt[Xt[:,bf]<=t,:]; Yt1=Yt[Xt[:,bf]<=t]
        Xv1=Xv[Xv[:,bf]<=t,:]; Yv1=Yv[Xv[:,bf]<=t]
        Xt2=Xt[Xt[:,bf]>t,:]; Yt2=Yt[Xt[:,bf]>t]
        Xv2=Xv[Xv[:,bf]>t,:]; Yv2=Yv[Xv[:,bf]>t]
    
        if (len(Yt1)==0)|(len(Yt2)==0):
            tree=LabelName[most(Yt)]
        else:
            tree[nt]['<=%.4f'%t]=CreatTree(Xt1,Yt1,Xv1,Yv1,FeatureName,LabelName,considerfeature.copy(),DataType,rule,cut,equalcut,deep+1,normal,lamuda,equalchoose,order)
            tree[nt]['>%.4f'%t]=CreatTree(Xt2,Yt2,Xv2,Yv2,FeatureName,LabelName,considerfeature.copy(),DataType,rule,cut,equalcut,deep+1,normal,lamuda,equalchoose,order)
    
    if cut=='after' and deep==0: #后剪枝
        tree,right=aftercut(tree,Xt,Yt,Xv,Yv,FeatureName,LabelName,DataType,equalcut)
    return tree